JARRIGE-DOMITILLE_evaluation-m6-2021

Introduction :

L’objectif de ce projet est de poursuivre les études que nous avons préalablement réalisées au cours du module 3 sur les données omiques issues de cet article :

Pavkovic, M., Pantano, L., Gerlach, C.V. et al. Multi omics analysis of fibrotic kidneys in two mouse models. Sci Data 6, 92 (2019). https://doi.org/10.1038/s41597-019-0095-5

Les auteurs étudient la fibrose rénale chez des souris. Ils utilisent deux modèles :

• un où la nephropathie est induite par un traitement chirurgical, irreversible

• un traitement à l’acide folique qui est reversible.

Les auteurs ont prélevé à différents temps après le traitement des échantillons de protéines et d’ARN de rein de souris.

Nous travaillerons sur le modèle FA (folic acid) en cherchant à mettre en relation les deux types de données omiques produites : transcriptomique et protéomique.

Je ne m’attarderai pas sur les analyses descriptives et exploratoires sur un seul type de données, car ces approches ont été déjà abordées lors des différents TP et le projet du module 3.

Note : il était difficile pour moi de me restreidre à deux figures pour certaines des analyses, puisque celles-ci servent souvent à guider mes choix et à les justifier. J’espère que vous ne m’en tiendrez pas rigueur…

setwd("/shared/projects/dubii2021/djarrige/m6-bioinfo-integr/mini-projet/")

mkdir -p data results

librairies = c("DESeq2",
               "dplyr",
               "VennDiagram",
               "tidyverse",
               "knitr",
               "WGCNA")

for (lib in librairies){
  if (!require(lib, character.only = TRUE)){
    install.packages(lib)
  }
  require(lib, character.only = TRUE)
}

required_Bioconductor <- c("mixOmics")

if (!requireNamespace("BiocManager", quietly = TRUE)){
  install.packages("BiocManager")
}

for (lib in required_Bioconductor) {
  if (!require(lib, character.only = TRUE)){
    BiocManager::install(lib)
  }
  require(lib, character.only = TRUE)
}

Préparation des données

# I download the raw data on the fa model and their corresponding metadata

url_base <- "https://github.com/DU-Bii/module-3-Stat-R/raw/master/stat-R_2021/data/pavkovic_2019/"
suffixes <- c("fa_raw_counts.tsv.gz",
              "fa_transcriptome_metadata.tsv",
              "pfa_model_counts.tsv.gz",
              "pfa_proteome_metadata.tsv")

for (suffix in suffixes){
  file_path <- paste0("data/", suffix)
  url <- paste0(url_base, suffix)
  if (!file.exists(file_path)){
    download.file(url, destfile = file_path)
  }
}

fa_raw <- read.csv("data/fa_raw_counts.tsv.gz", sep="\t", row.names=1)
pfa_raw <- read.csv("data/pfa_model_counts.tsv.gz", sep="\t")

fa_meta <- read.csv("data/fa_transcriptome_metadata.tsv", sep="\t", row.names=1)
pfa_meta <- read.csv("data/pfa_proteome_metadata.tsv", sep="\t", row.names=1)

Dans les données de protéomique, des id sont dupliqués. Pour pouvoir les utiliser comme row.names il va falloir les renommer.

id <- pfa_raw$id  # extracts the id 
id <- make.unique(id) # make them unique by appending numbers after duplicates
row.names(pfa_raw) <- id # use those unique id as rownames for the proteomics data
pfa_raw <- pfa_raw[, -1] # I take out the former id column

Les échantillons ne sont pas dans le même ordre dans les metadonnées transcriptomiques et les données. Je corrige cela !

ref_order <- paste0(rep(c("normal", "day1", "day2", "day3", "day7", "day14"), each=3), "_", 1:3)

fa_meta <- fa_meta[match(ref_order, fa_meta$sampleName),]
fa_raw <- fa_raw[, match(ref_order, colnames(fa_raw))]

Analyse d’expression différentielle :

Consignes :

Analyse d’expression différentielle pour les données de protéomique et transcriptomique => identifier les gènes/protéines significativement différentiellement exprimés dans le modèle FA en comparant Day 7 à Day 0.

Je vais utiliser DESeq2 pour traiter les données transcriptomiques. Avant le traitement par DESeq2 il faut :

-que les données de comptage soient brutes et entières (je vais considérer ici que les données dites raw des auteurs le sont bien réellement…)

-transformer le dataframe de comptage en matrice

Traitement des données de transcriptomique

fa_raw_matrix <- as.matrix(fa_raw)
fa_raw_matrix <- round(fa_raw_matrix)

fa_meta_simplified <- fa_meta[, c("sampleName", "condition")]

fa_deseq <- DESeqDataSetFromMatrix(countData = fa_raw_matrix, 
                                   colData = fa_meta,
                                   design = ~condition)

## converting counts to integer mode

Retrait des gènes pas ou quasi pas exprimés dans les données transcriptomiques, on ne garde que ceux dont le total de comptages est supérieur à 10.

# taking out low expression genes
fa_deseq <- fa_deseq[ rowSums(counts(fa_deseq)) > 10, ]

La fonction DESeq() réalise plusieurs étapes :

• l’estimation des size factors, pour normaliser et prendre en compte la taille des libraries dans l’analyse

• l’estimation de la dispersion pour chaque gène

• un fit sur un modèle linéaire

Ensuite on s’interessera qu’aux conditions “normal” et “day7” en les extrayant des résultats.

Analyse différentielle DESeq2 transcriptomique

fa_deseq <- DESeq(fa_deseq)

# Extracting differential analyse results for conditions day7 on normal
res_DESeq <- results(fa_deseq, contrast = c("condition", "day7", "normal"))

# selecting genes differentially expressed with an adjusted (Benjamini-Hochberg) p value under 5%
significant_genes <- data.frame(id = row.names(res_DESeq[which(res_DESeq$padj < 0.05),]),
                                log2FoldChange_transcriptomics = res_DESeq[which(res_DESeq$padj < 0.05), "log2FoldChange"])

row.names(significant_genes) <- significant_genes$id

J’essaie DESeq2 également avec les données de protéomique, car celui-ci peut en théorie les traiter (une mention discrète dans la documentation) cependant il n’est pas construit pour les données de protéomique ce qui pourrait potentiellement biaiser le résultat. Nous allons essayer de voir cela.

Traitement des données de protéomique

pfa_raw_matrix <- as.matrix(pfa_raw)
pfa_raw_matrix <- round(pfa_raw_matrix)

pfa_deseq <- DESeqDataSetFromMatrix(countData = pfa_raw_matrix,
                                    colData = pfa_meta,
                                    design = ~ condition)

## converting counts to integer mode

Analyse différentielle DESeq2 protéomique :

pfa_deseq <- DESeq(pfa_deseq)
res_prot_DESeq2 <- results(pfa_deseq, contrast = c("condition", "day7", "normal"))

# Saving file for network analysis
tmp_df <- as.data.frame(res_prot_DESeq2)
tmp_df <- add_column(.data = tmp_df, gene_id = row.names(tmp_df), .before = 1)

write_csv(tmp_df, file = "results/protein_logFoldChange_7-0.csv")
remove(tmp_df)

significant_proteins <- data.frame(id = row.names(res_prot_DESeq2[which(res_prot_DESeq2$padj < 0.05), ]),
                                   log2FoldChange_proteomics = res_prot_DESeq2[which(res_prot_DESeq2$padj < 0.05), "log2FoldChange"])
row.names(significant_proteins) <-significant_proteins$id

Visualisation des résultats :

Les gènes ou protéines significativement (p value adjusted < 0.05%) différentiellement exprimés sont colorés en bleu.

par(mfrow=c(1,2))
plotMA(res_DESeq, las=1, cex.main = 1, ylim=c(-8,8),
       main = "Differentially expressed genes")
plotMA(res_prot_DESeq2, las=1, cex.main = 1, ylim=c(-8,8),
       main = "Differentially expressed proteins")

Differentially expressed genes and proteins.

significant_all <- full_join(significant_genes, significant_proteins, by="id")
row.names(significant_all) <- significant_all$id

print(paste0("Differentially expressed RNAs: ", dim(significant_genes)[1]))

## [1] "Differentially expressed RNAs: 1996"

print(paste0("Differentially expressed proteins: ", dim(significant_proteins)[1]))

## [1] "Differentially expressed proteins: 4462"

print(paste0("Genes with both differentially expressed RNA and protein: ", dim(na.omit(significant_all))[1]))

## [1] "Genes with both differentially expressed RNA and protein: 841"

kable(na.omit(significant_all)[1:10,], caption = "Some examples:")

Some examples:

	id	log2FoldChange_transcriptomics	log2FoldChange_proteomics
ENSMUSG00000000078	ENSMUSG00000000078	2.6721831	-1.0822107
ENSMUSG00000000154	ENSMUSG00000000154	-1.8348021	-2.3124483
ENSMUSG00000000171	ENSMUSG00000000171	-1.0890533	-1.0857821
ENSMUSG00000000563	ENSMUSG00000000563	-1.3407229	-1.1570848
ENSMUSG00000000673	ENSMUSG00000000673	-2.3409525	-2.1790929
ENSMUSG00000000876	ENSMUSG00000000876	-0.8951156	-1.4648656
ENSMUSG00000000915	ENSMUSG00000000915	1.6325695	0.2170728
ENSMUSG00000001025	ENSMUSG00000001025	2.4919169	1.6601354
ENSMUSG00000001089	ENSMUSG00000001089	2.3324185	0.5589445
ENSMUSG00000001100	ENSMUSG00000001100	-1.5391918	-0.8996146

strange_genes <-  c(na.omit(significant_all[which(significant_all$log2FoldChange_transcriptomics > 0 & significant_all$log2FoldChange_proteomics < 0), "id"]), na.omit(significant_all[which(significant_all$log2FoldChange_transcriptomics < 0 & significant_all$log2FoldChange_proteomics > 0), "id"]))


kable(significant_all[strange_genes, 2:3], caption = "Unconsistancy between RNA and protein")

Unconsistancy between RNA and protein

	log2FoldChange_transcriptomics	log2FoldChange_proteomics
ENSMUSG00000000078	2.6721831	-1.0822107
ENSMUSG00000001247	1.4046401	-0.2860428
ENSMUSG00000020744	2.7882125	-0.9110466
ENSMUSG00000024096	1.6492218	-0.3060593
ENSMUSG00000025742	1.1677440	-0.7607891
ENSMUSG00000040264	2.1397001	-0.7939093
ENSMUSG00000041797	1.5940262	-0.2789294
ENSMUSG00000002103	-1.7919769	0.8207154
ENSMUSG00000017210	-2.4336964	0.3247741
ENSMUSG00000020064	-1.0303276	0.1909243
ENSMUSG00000020180	-1.8407726	0.4060388
ENSMUSG00000021468	-1.4980361	0.2506465
ENSMUSG00000021660	-0.8394688	0.2133396
ENSMUSG00000021975	-1.4776663	0.4075119
ENSMUSG00000022205	-1.4609761	0.7303203
ENSMUSG00000028189	-1.4827506	0.4399862
ENSMUSG00000031660	-2.6095592	0.6029440
ENSMUSG00000031848	-1.8760340	0.5725588
ENSMUSG00000036181	-1.4033241	0.7846991
ENSMUSG00000042626	-1.9752141	0.4720450
ENSMUSG00000045996	-1.3126510	0.8334348
ENSMUSG00000051695	-1.1590639	0.3449563
ENSMUSG00000059796	-1.7387242	0.6464077
ENSMUSG00000060639	-1.7564435	0.7010975
ENSMUSG00000061048	-2.8027373	2.2505599

print(paste0("Unconsistant genes: ", (length(strange_genes) / nrow(significant_all) * 100), " %"))

## [1] "Unconsistant genes: 0.445077443475165 %"

Le diagramme de Venn ci-dessous présente l’intersection entre les gènes différentiellement exprimés au niveau ARN et ceux au niveau protéine, entre la condition normale et day7. Environ 840 gènes voient leur niveau significativement différent pour les deux types de molécule.

On pourrait ensuite par exemple réaliser une analyse d’enrichissement fonctionnel sur cette liste.

Attention cependant, comme on peut le constater dans les exemples ci-dessus quelques gènes peuvent être surexprimés en ARN et sous exprimés en protéines (ou inversement) ! Heureusement ces gènes aux expressions contradictoires sont rares ! Difficile de dire sans informations complémentaires s’il pourrait s’agir d’une réalité biologique (par exemple la protéine pourrait être soumise à une dégradation très intense et le rein pourrait tenter de compenser en produisant beaucoup d’ARNm), d’un biais technique, ou bien si l’utilisation de DESeq2 sur des données de protéomique pourrait entrainer des artéfacts… Beaucoup de protéines sont jugées différentiellement exprimées par DESeq2 en comparaison aux ARN.

area1 <- dim(significant_proteins)[1]
area2 <- dim(significant_genes)[1]
cross <- dim(na.omit(significant_all))[1]

venn.plot <- draw.pairwise.venn(area1 = area1, 
                                area2 = area2, 
                                cross.area = cross,
                                fontfamily = "sans",
                                category = c("significant proteomics", "significant transcriptomics"),
                                fill = c(alpha("magenta",0.3), alpha('cyan',0.3)),
                                lwd=0,
                                cat.col = c("magenta", "cyan"),
                                cat.pos = c(0, 0),
                                cat.dist = c(0.02, 0.087),
                                cat.fontfamily = "sans")
grid.draw(venn.plot)

Differentially expressed genes and proteins, overlap.

grid.newpage()

Analyse multi-omique

Consignes :

Analyse multi-omique (transcripto + protéo) avec, au choix, MOFA, mixOmics, mixKernel ou d’autres outils de factorisation multi-matrices. Vous pouvez soit vous focaliser sur un time point, soit intégrer les différents time points, en partant des données normalisées fournies dans le matériel supplémentaire du papier.

J’utilise mixOmics. Cependant je vais essayer de normaliser les données moi-même (en m’inspirant de celles utilisées dans le module 3) et de regarder tous les timempoints possibles (à part le jour 3 qui est absent des données de protéomiques).

Avant d’explorer nos données avec mixOmics il faut :

• les normaliser

• les filtrer (pour avoir maximum environ 10 000 gènes ou protéines)

• passer les échantillons en lignes et les gènes/protéines en colonnes

• avoir les mêmes échantillons dans les deux analyses

Pour les données de protéomique, on a un nombre d’observables raisonnable (8044), je ne filtrerai donc pas de protéines.

En revanche pour les données transcriptomiques, on a beaucoup d’ARNs (46679), il vaut mieux les filtrer !

Traitement des données

Retrait des gènes pas ou quasi pas exprimés dans les données transcriptomiques, on ne garde que ceux dont la somme des comptage est supérieure à 10.

fa_filtered <- fa_raw[which(rowSums(fa_raw) > 10), ]

Normalisation inter-échantillon sur les troisièmes quartiles

third_quantiles <- apply(fa_filtered, 2 , quantile, p=0.75, na.rm = TRUE)
global_third_quantile <- quantile(unlist(fa_filtered), p=0.75)
size_factors <- 1 / third_quantiles * global_third_quantile


fa_norm <- t(t(fa_filtered) * size_factors)

Transformation en log2

fa_log2 <- log2(fa_norm + 1)

#boxplot(fa_log2, col=fa_meta$color)

Calculs de statistiques sur les gènes pour ne garder que les gènes les plus exprimés dans l’analyse.

gene_stats <- data.frame(row.names = row.names(fa_log2))
gene_stats$mean <- apply(fa_log2, 1, mean)
gene_stats$var <- apply(fa_log2, 1, var)

Selection des gènes exprimés les plus variables.

gene_kept <- gene_stats[which((gene_stats$mean > 4) & (gene_stats$var > 1.5)), ]
fa_log2 <- fa_log2[row.names(gene_kept),]

#boxplot(fa_log2, col=fa_meta$color)

Normalisation inter-échantillons des données de protéomique

third_quantiles <- apply(pfa_raw, 2 , quantile, p=0.75, na.rm = TRUE)
global_third_quantile <- quantile(unlist(pfa_raw), p=0.75)
size_factors <- 1 / third_quantiles * global_third_quantile


pfa_norm <- t(t(pfa_raw) * size_factors)

Transformation des données de protéomique en log2

pfa_log2 <- log2(pfa_norm + 1)

#boxplot(pfa_log2, col=pfa_meta$color)

Je vais retirer les réplicats biologique n°3, afin d’avoir les même individus dans les deux jeux de données. Je retire également les timepoints du jour 3, qui sont absents du jeu protéomique.

fa_mix <- fa_log2[, -c(3, 6, 9, 10, 11, 12, 15, 18)]
meta_mix <- fa_meta[-c(3, 6, 9, 10, 11, 12, 15, 18),-1]

multi_dataset <- list(genes = t(fa_mix),
                      proteins = t(pfa_log2),
                      conditions = meta_mix)

Sparse Partial to Least Squares (sPLS)

Cette méthode est non-supervisée, les échantillons ne sont pas séparés en groupes. La sparse PLS va intégrer les données de transcriptomique et de protéomique tout en sélectionnant les variables (gènes ou protéines) les plus covariants ou anti-covariants dans les deux jeux de données omiques.

mix_sPLS <- spls(X = multi_dataset$genes,
                 Y = multi_dataset$proteins,
                 keepX= c(25, 25),
                 keepY= c(25, 25))

plotIndiv(mix_sPLS, rep.space = 'XY-variate', cex=3.5,
          group = multi_dataset$conditions$condition,
          title = "Individuals Sparse PLS")

sPLS on transcriptomic and proteomic datasets

Ici, je demande à la méthode de garder 25 variables de chaque jeu de données par composante. Je demande d’inclure deux composantes dans le modèle.

Le plot présente une synthèse de l’analyse des deux jeux de données.

On constate que la sPLS nous permet non seulement d’étudier la relation entre nos deux jeux de données, mais ici permet de séparer “naturellement” (sans supervision) les échantillons en “clusters” distincts. Les normaux isolés, les conditions peu après l’injection proches et celles plus tardives ensemble.

par(mfrow=c(1,2))
plotLoadings(mix_sPLS, size.title=1.2, comp=1, 
             title = "Loadings on comp 1",
             subtitle = c("RNAs", "Proteins"),
             size.subtitle = 1)

sPLS: Loadings on genes and proteins.

plotLoadings(mix_sPLS, size.title=1.2, comp=2,
             title = "Loadings on comp 2",
             subtitle = c("RNAs", "Proteins"),
             size.subtitle = 1)

sPLS: Loadings on genes and proteins.

Ici les loading représentent quels gènes contribuent le plus à la covariance entre les deux jeux de données omiques. Le bloc X représente les gènes dont l’expression caractérise la relation entre données transcriptomiques et protéomiques, le bloc Y les protéines.

On obtient donc des listes de gènes et protéines qui pourront caractériser les échantillons.

Réseau de co-expression avec WGCNA

Consignes :

Construction de réseau avec WGCNA à partir des données transcriptomiques.

Reconstruct the co-expression network from all the time points of the FA transcriptomics data. Propose to filter and remove all the zero expressed genes, the NAs and the less informative genes from the transcriptomics data. (I remove all the genes that are not expressed in at least 9 out of the 18 conditions (expression > 1 TPM in 9) and then filter with the coefficient of variation > 0.75).

Then apply the first part of the network reconstruction steps as we saw them on the WGCNA course until the module predictions.

Instead of using WGCNA’s module prediction routines, apply a universal threshold of 0.5 on the adjacency matrix, and obtain an adjacency matrix that is reduced in size. This is the network. Import it to Cytoscape with aMatReader plugin. Visualize, analyze the network and superimpose the proteomics data on it.

Filtrage des données de transcriptomique

# Computing gene statistics
stats_wgcna <- data.frame(mean = apply(fa_raw, 1, mean),
                          sd = apply(fa_raw, 1, sd),
                          row.names = row.names(fa_raw))

stats_wgcna$coef_var <- (stats_wgcna$sd / stats_wgcna$mean)
stats_wgcna$null = apply(fa_raw == 0, 1, sum, na.rm = TRUE)


# filtering at least expressed in 9 samples and coefficient of variation superior to 0.75

fa_wgcna <- fa_raw[ row.names(stats_wgcna[ which(stats_wgcna$null < 9 & stats_wgcna$coef_var > 0.75), ]), ]

Je réalise aussi une normalisation et une transformation log2, car sinon le clustering est très différent de celui obtenu lors du mini-projet du module 3.

third_quantiles <- apply(fa_wgcna, 2 , quantile, p=0.75, na.rm = TRUE)
global_third_quantile <- quantile(unlist(fa_wgcna), p=0.75)
size_factors <- 1 / third_quantiles * global_third_quantile

fa_wgcna <- t(t(fa_wgcna) * size_factors)

fa_wgcna <- log2(fa_wgcna + 1)

sample_tree = hclust(dist(t(fa_wgcna)), method = "complete")

plot(sample_tree, main = "Sample clustering to detect outliers", 
     xlab="", cex.main = 1)

abline(h = 250, col = "red")

# For me day2_2 is an outlier, I take it out.
fa_wgcna <- fa_wgcna[,-8]

sample_tree2 <- hclust(dist(t(fa_wgcna)), method = "complete")

Je retire l’échantillon day2_2 qui est un outlier

Maintenant je crée un data.frame avec des données binaires, reflétant le temps de prélèvement.

traits_data <- data.frame(names = fa_meta$sampleName)

traits_data$normal <- as.integer(startsWith(traits_data$names, "normal"))
traits_data$day1 <- as.integer(startsWith(traits_data$names, "day1_"))
traits_data$day2 <- as.integer(startsWith(traits_data$names, "day2"))
traits_data$day3 <- as.integer(startsWith(traits_data$names, "day3"))
traits_data$day7 <- as.integer(startsWith(traits_data$names, "day7"))
traits_data$day14 <- as.integer(startsWith(traits_data$names, "day14"))

row.names(traits_data) <- traits_data$names
traits_data <- traits_data[,-1]
traits_data <- traits_data[-8,]

Traits et échantillons

Les différents temps de prélèvement corrèlent plutôt bien avec ce qui est attendu pour les échantillons.

# Convert traits to a color representation: red means yes
traits_colors = numbers2colors(traits_data, signed = FALSE);

# Plot the sample dendrogram and the colors underneath.
plotDendroAndColors(sample_tree2, traits_colors,
                    groupLabels = names(traits_data),
                    main = "Sample dendrogram and traits heatmap")

Je transpose la matrice de données d’expression pour la suite de l’analyse.

fa_wgcna <- t(fa_wgcna)

WGCNA construction du réseau

Recherche d’une puissance adaptée

save(fa_wgcna, traits_data, file = "results/fa_transcriptomics_wgcna.RData")

allowWGCNAThreads()

## Allowing multi-threading with up to 56 threads.

# Load the data saved in the first part
lnames = load(file = "results/fa_transcriptomics_wgcna.RData");


# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 12, to = 20, by = 2))

# Call the network topology analysis function
sft = pickSoftThreshold(fa_wgcna, powerVector = powers, verbose = 0)

##    Power SFT.R.sq slope truncated.R.sq  mean.k. median.k. max.k.
## 1      1    0.180 -3.44          0.683 2580.000  2530.000 3480.0
## 2      2    0.765 -4.05          0.900  849.000   799.000 1560.0
## 3      3    0.924 -3.68          0.965  346.000   310.000  862.0
## 4      4    0.955 -3.34          0.971  163.000   138.000  547.0
## 5      5    0.956 -3.07          0.967   84.600    68.000  380.0
## 6      6    0.963 -2.81          0.969   47.600    35.900  281.0
## 7      7    0.965 -2.62          0.972   28.600    20.100  217.0
## 8      8    0.956 -2.48          0.963   18.000    11.800  174.0
## 9      9    0.952 -2.37          0.962   11.900     7.220  143.0
## 10    10    0.957 -2.24          0.968    8.150     4.560  120.0
## 11    12    0.960 -2.05          0.976    4.190     1.980   87.8
## 12    14    0.945 -1.95          0.969    2.380     0.939   67.4
## 13    16    0.933 -1.86          0.969    1.460     0.478   53.4
## 14    18    0.934 -1.78          0.974    0.954     0.258   43.3
## 15    20    0.933 -1.70          0.972    0.655     0.146   35.7

# Plot the results:
par(mfrow = c(1, 2));
options(repr.plot.width = 14, repr.plot.height = 10);

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit,signed R^2", type = "n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels = powers, cex = 0.9, col = "red");
# this line corresponds to using an R^2 cut-off of h
abline(h = 0.90, col = "red")

# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, cex = 0.9, col = "red")

Soft thresholding

Je prends le power 6 pour la suite du pipeline.

Prédiction de modules

net = blockwiseModules(fa_wgcna, power = 6,
                       TOMType = "signed", minModuleSize = 30,
                       reassignThreshold = 1e-6, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE, nThreads = 8,
                       saveTOMFileBase = "results/fa_transcriptomic_TOM",
                       verbose = 0)

# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
#mergedColors
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

Dendrogram and modules.

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
     file = "results/fa_transcriptomic_networkConstruction-auto.RData")



lnames = load(file = "results/fa_transcriptomics_wgcna.RData");

# Load network data saved in the second part.
lnames = load(file = "results/fa_transcriptomic_networkConstruction-auto.RData");

# Define numbers of genes and samples
nGenes = ncol(fa_wgcna);
nSamples = nrow(fa_wgcna);

# Recalculate MEs with color labels
MEs0 = moduleEigengenes(fa_wgcna, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, traits_data, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

options(repr.plot.width=16, repr.plot.height=12)
# Will display correlations and their p-values
textMatrix =  paste(signif(moduleTraitCor, 2), "\n(",
                           signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 11, 1, 0));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(traits_data),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               cex.legendLabel = 0.6,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))

Module-trait relationships

J’ai tenté de comparer l’utilisation du TOM sans prédiction de module, représentant l’ensemble des gènes, avec le modTOM qui ne contient qu’une sélection de deux modules d’intérêt.

J’ai rencontré des difficultés avec l’utilisation du threshold, mes données ne dépassant jamais le niveau demandé de 0.5 (voir histogrammes) malgré l’essai de nombreux paramètres, puissances et méthodes.

J’ai alors supposé que ce threshold était basé sur l’adjacency et pas directement la corrélation, et en ai calculé un en suivant la formule de calcul de la méthode “signed” sur une corrélation de 0.5 et un power de 6, soit environ 0.18.

# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(fa_wgcna, power = 6,
                            networkType = "signed",
                            verbose = 0)

# Select modules
# I pick one module strongly expressed in normal samples and one in day 7 samples
modules = c("yellow", "white");  

# Select module probes
probes = colnames(fa_wgcna)
inModule = is.finite(match(moduleColors, modules))
modProbes = probes[inModule]

# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule];
dimnames(modTOM) = list(modProbes, modProbes)

####
#Adjacency values visualization to decide on threshold
par(mfrow=c(1,2))
hist(unlist(TOM), breaks = 100, cex.main = 0.8)
hist(unlist(modTOM), breaks = 100, cex.main = 0.8)

# threshold : "adjacency threshold for including edges in the output." (from the doc)
# With "signed" type adjacency -> adjacency = (0.5 * (1+cor) )^power
# if correlation = 0.5 -> we get adjacency = 0.18

Exportation des données pour Cytoscape

J’essaye de produire un sous réseau avec la sélection de deux modules puis le réseau complet.

# Export the network into edge and node list files Cytoscape can read
### Here I keep the whole network
cyt = exportNetworkToCytoscape(TOM,
  edgeFile = paste("results/CytoscapeInput-edges", ".txt", sep=""),
  nodeFile = paste("results/CytoscapeInput-nodes", ".txt", sep=""),
  weighted = TRUE,
  threshold = 0.18, # I modified the threshold here
  nodeNames = probes,
  #nodeAttr = moduleColors[inModule]
  )

# Export the network into edge and node list files Cytoscape can read
### Here I export only two gene modules of interest
cyt = exportNetworkToCytoscape(modTOM,
  edgeFile = paste("results/CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
  nodeFile = paste("results/CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
  weighted = TRUE,
  threshold = 0.18, # I modified the threshold here
  nodeNames = modProbes,
  nodeAttr = moduleColors[inModule]
  )

Visualisation du réseau dans Cytoscape

Consigne :

Colorez dans le réseau choisi les noeuds en fonction des données de protéomiques avec un gradient de couleur correspondant au fold-change des données de protéomique.

Sous-réseau (2 modules)

Visualisons d’abord le sous-réseau avec les deux modules sélectionnés.

Les gènes se séparent en deux réseaux. J’ai coloré les nodes en fonction de leur appartenance au module yellow (plus exprimé en condition normale) ou white (plus exprimé en condition day7). Etrangement je n’observe pas de séparation évidente des gènes en fonction de leur appartenance aux modules.

reseau_2_modules

Réseau complet

Maintenant, visualisons le réseau complet des gènes variables :

Dans cette version brute on retrouve la séparation en deux réseaux isolés. De plus on constate que celui de gauche semble avoir deux parties moins fortement correlées.

reseau_entier

Puis, je n’ai conservé que les gènes pour lesquels une information protéique était disponible. Ensuite, j’ai coloré les nodes en fonction du log fold change entre les conditions normales et jour7.

On constate ici une séparation franche des gènes plus exprimés en condition day7 qu’en condition normal (en rouge) de ceux moins exprimés qu’en condition normale (en bleu).

Le réseau de droite regroupe des gènes plus exprimés lors de la nephropathie. Celui de gauche est intéressant ; il contient un “pôle” de gènes moins exprimés en condition day7 qu’en normal (reflétant certainement l’activité normale d’un rein sain) et un petit “pôle” de gènes à l’inverse plus exprimés.

Nos données protéomiques donnent une cohérence supplémentaire à ce réseau inféré avec WGCNA. Aussi, je pense que l’utilisation de DESeq2 pour les données de spectrométrie de masse n’est peut-être pas idéale, mais donne un résulat final globalement cohérent avec les données de RNA-seq.

reseau_entier_prot

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /shared/ifbstor1/software/miniconda/envs/r-4.0.3/lib/libopenblasp-r0.3.10.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] mixOmics_6.14.0             lattice_0.20-41            
##  [3] MASS_7.3-53.1               WGCNA_1.70-3               
##  [5] fastcluster_1.1.25          dynamicTreeCut_1.63-1      
##  [7] knitr_1.33                  forcats_0.5.1              
##  [9] stringr_1.4.0               purrr_0.3.4                
## [11] readr_1.4.0                 tidyr_1.1.3                
## [13] tibble_3.1.2                ggplot2_3.3.3              
## [15] tidyverse_1.3.0             VennDiagram_1.6.20         
## [17] futile.logger_1.4.3         dplyr_1.0.6                
## [19] DESeq2_1.30.1               SummarizedExperiment_1.20.0
## [21] Biobase_2.50.0              MatrixGenerics_1.2.1       
## [23] matrixStats_0.58.0          GenomicRanges_1.42.0       
## [25] GenomeInfoDb_1.26.7         IRanges_2.24.1             
## [27] S4Vectors_0.28.1            BiocGenerics_0.36.1        
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1           backports_1.2.1        Hmisc_4.5-0           
##   [4] plyr_1.8.6             igraph_1.2.6           splines_4.0.3         
##   [7] BiocParallel_1.24.1    digest_0.6.27          foreach_1.5.1         
##  [10] htmltools_0.5.1.1      GO.db_3.12.1           fansi_0.4.2           
##  [13] magrittr_2.0.1         checkmate_2.0.0        memoise_2.0.0         
##  [16] cluster_2.1.1          doParallel_1.0.16      annotate_1.68.0       
##  [19] modelr_0.1.8           rARPACK_0.11-0         jpeg_0.1-8.1          
##  [22] colorspace_2.0-1       blob_1.2.1             rvest_1.0.0           
##  [25] ggrepel_0.9.1          haven_2.3.1            xfun_0.23             
##  [28] crayon_1.4.1           RCurl_1.98-1.3         jsonlite_1.7.2        
##  [31] genefilter_1.72.1      impute_1.64.0          survival_3.2-10       
##  [34] iterators_1.0.13       glue_1.4.2             gtable_0.3.0          
##  [37] zlibbioc_1.36.0        XVector_0.30.0         DelayedArray_0.16.3   
##  [40] scales_1.1.1           futile.options_1.0.1   DBI_1.1.1             
##  [43] Rcpp_1.0.6             xtable_1.8-4           htmlTable_2.2.1       
##  [46] foreign_0.8-81         bit_4.0.4              preprocessCore_1.52.1 
##  [49] Formula_1.2-4          htmlwidgets_1.5.3      httr_1.4.2            
##  [52] RColorBrewer_1.1-2     ellipsis_0.3.2         farver_2.1.0          
##  [55] pkgconfig_2.0.3        XML_3.99-0.6           nnet_7.3-15           
##  [58] sass_0.4.0             dbplyr_2.1.1           locfit_1.5-9.4        
##  [61] utf8_1.2.1             labeling_0.4.2         tidyselect_1.1.1      
##  [64] rlang_0.4.11           reshape2_1.4.4         AnnotationDbi_1.52.0  
##  [67] munsell_0.5.0          cellranger_1.1.0       tools_4.0.3           
##  [70] cachem_1.0.5           cli_2.5.0              generics_0.1.0        
##  [73] RSQLite_2.2.7          broom_0.7.5            evaluate_0.14         
##  [76] fastmap_1.1.0          yaml_2.2.1             bit64_4.0.5           
##  [79] fs_1.5.0               formatR_1.9            xml2_1.3.2            
##  [82] compiler_4.0.3         rstudioapi_0.13        png_0.1-7             
##  [85] reprex_1.0.0           geneplotter_1.68.0     bslib_0.2.5.1         
##  [88] stringi_1.6.2          highr_0.9              RSpectra_0.16-0       
##  [91] Matrix_1.3-2           vctrs_0.3.8            pillar_1.6.1          
##  [94] lifecycle_1.0.0        BiocManager_1.30.10    jquerylib_0.1.4       
##  [97] data.table_1.14.0      bitops_1.0-7           corpcor_1.6.9         
## [100] R6_2.5.0               latticeExtra_0.6-29    gridExtra_2.3         
## [103] codetools_0.2-18       lambda.r_1.2.4         assertthat_0.2.1      
## [106] withr_2.4.2            GenomeInfoDbData_1.2.4 hms_1.1.0             
## [109] rpart_4.1-15           rmarkdown_2.8          lubridate_1.7.10      
## [112] base64enc_0.1-3        ellipse_0.4.2